Introduction à l’analyse des données
1 Markdown 101
Markdown est un langage de balisage et peut être utilisé pour donner aux documents texte tels que R Markdown une structure comme word ou latex. La syntaxe est assez simple.
quelques exemples (tiré de https://www.markdownguide.org/basic-syntax):
# Heading level 1
Titre de section
## Heading level 2
Titre de sous-section
### Heading level 3
Titre de sous-sous section
etc..
__un texte en gras__
_un texte en italique_
~~un text barré~~
ie
un texte en gras
un texte en italique
un text barré
Pour plus d’informations, consultez par exemple le lien de rstudio.
1.1 LaTeX
R Markdown est encore plus puissant, on peut utilisier \(\LaTeX\), soit inline comme \(y_i=\beta_0 + x_i\beta_1\) ou comme equation: \[ \hat{\beta}=(\boldsymbol{X}^\top\boldsymbol{X})^{-1}\boldsymbol{X}^\top\boldsymbol{y} \] vous pouvez même inclure des fichiers .tex entiers.
1.2 Figures
Imuages (ce n’est pas un plot) peuvent également être facilement incorporés avec 
1.3 Knitting
Pour les documents de base, le yaml est ce qui est important:
rmdformats::robobook c’est seulement pour rendre le ficher légèrement plus joli
1.4 Pourquoi?
- Projets pour votre thèse de baccalauréat (Peut être entièrement reproductible)
- après votre baccalauréat
- Industrie-1: Vous ferez beaucoup de rapports, souvent répétitifs, ce qui n’est pas très fun.. mais vous pouvez automatiser la plupart avec markdown
- Industrie-2: Particulierement en tant que junior, vous aurez besoin de faire des recherches, Rmarkdown est parfait pour les présenter
- Université: Si vous faites une maitrise - Vous devrez probablement coder et remettre votre code
- Vous pouvez également utiliser python à partir des chunks de rmarkdown
2 Manipulation des Données
Si vous travaillez avec des données, vous passerez beaucoup de temps à préparer vos données avant de modeliser (ensembles de données “propres” n’existent généralement pas dans le monde réel). Dans le cadre de ce cours, nous n’examinerons que les principes de base des “tidy data”.
L’idée générale est de mettre de l’ordre dans nos données de façon à ce que nous ayons :
- une observation par ligne
- Une information par cellule
- Un type d’information par colonne
Pour plus il y a le résumé du papier de Hadley Wickham, vous pouvez consulter ce lien
Dans ce qui suit, nous envisagerons toujours plusieurs façons de résoudre un problème. Vous êtes libre de choisir celle que vous correspond le plus, bien que je recommande personnellement l’écosystème “tidyverse”.
Pour des raisons esthétiques, il est souvent recommandé d’utiliser l’opérateur pipe (%>%). Sous R-Studio avec Cmd-Shift-M ou Ctrl-Shift-M. Ce qu’il fait fondamentalement est de prendre la ligne précédente comme premier argument dans la ligne courante.
# meme chose
my_numbers <- c(0,1,2,3)
print(sum(my_numbers))## [1] 6
print(my_numbers %>% sum())## [1] 6
cette logique simplifiera grandement vos flux de travail et rendra le code plus lisible.
2.1 SQL-Like
SQL (Structured Query Language) - votre meilleur ami si vous travaillez avec des données dans une entreprise, mais également utile à plus petite échelle (aka. pour vos projets).
L’exercice est basé sur l’article : “Where is the Land of Opportunity ? The Geography of Intergenerational Mobility in the United States” (Raj Chetty, Nathaniel Hendren, Patrick Kline et Emmanuel Saez ; Quarterly Journal of Economics 129(4) : 1553-1623, 2014).Il fait partie d’un plus grand projet disponible à l’adresse suivante : http://www.equality-of-opportunity.org/.
Nous commençons par charger quelques données échantillons. Notez que “kable” dans la dernière ligne rend seulement les graphiques jolis.
library(haven) # si le format est .dta
library(readxl) # se le format est .xls / .xlsx
mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")
mobility_data %>%
head() %>%
kable()| index | cz | czname | e_rank_b | cs_race_bla | cs_race_theil_2000 | cs00_seg_inc | cs00_seg_inc_pov25 | cs00_seg_inc_aff75 | frac_traveltime_lt15 | hhinc00 | gini | inc_share_1perc | gini99 | frac_middleclass | taxrate | subcty_total_expenditure_pc | tax_st_diff_top20 | eitc_exposure | ccd_exp_tot | ccd_pup_tch_ratio | score_r | dropout_r | num_inst_pc | tuition | gradrate_r | cs_labforce | cs_elf_ind_man | d_tradeusch_pw_1990 | frac_worked1416 | mig_inflow | mig_outflow | cs_born_foreign | scap_ski90pcm | rel_tot | crime_violent | cs_fam_wkidsinglemom | cs_divorced | cs_married |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 100 | Johnson City | 43.58841 | 0.0693926 | 0.0576880 | 0.0408645 | -0.0171138 | 0.0266398 | 0.3844612 | 38209.54 | 0.4193436 | 5.695437 | 0.2528637 | 0.5840611 | 0.0158913 | 1404.822 | 1.0794937 | 0.0000780 | 6.577284 | 17.96812 | NA | NA | 0.0554321 | 1032.0113 | -0.2210098 | 0.6529477 | 0.0600338 | NA | 0.0148774 | 0.0417411 | 0.0559070 | 0.1604878 | NA | 0.5358987 | 0.0482720 | 0.2421846 | 0.1392865 | 0.5491035 |
| 1 | 200 | Morristown | 39.64703 | 0.4453178 | 0.0749721 | -0.0258754 | 0.0352996 | 0.0089049 | 0.3543061 | 25123.34 | 0.4867183 | 11.816536 | 0.3926104 | 0.3559214 | 0.0257373 | 1568.710 | 0.0198662 | -0.0128659 | 4.405095 | 17.79943 | -16.015344 | 0.0561887 | 0.0379247 | 970.0455 | -0.0416078 | 0.5273348 | 0.2086723 | 0.7670213 | 0.0223866 | 0.0181590 | 0.0137850 | 0.0232165 | -0.6127989 | 0.6469429 | 0.0316044 | 0.3499764 | 0.1039057 | 0.5426201 |
| 2 | 301 | Middlesborough | 46.93892 | 0.0151223 | 0.1100580 | 0.0954185 | 0.0139200 | 0.0527287 | 0.4793166 | 32701.43 | 0.3071970 | 11.810008 | 0.2077348 | 0.6675747 | 0.0116976 | 3115.946 | -0.0046168 | 0.0343010 | 6.609784 | 17.85102 | -9.137633 | 0.0246227 | 0.0634999 | 1395.0190 | -0.1716805 | 0.7066224 | 0.0362908 | 0.0274560 | 0.0274638 | 0.0309692 | 0.0421808 | 0.1318125 | -2.0635338 | 0.4222386 | 0.0511852 | 0.1831544 | 0.1402021 | 0.6224286 |
| 3 | 302 | Knoxville | NA | -0.0004112 | -0.0044561 | 0.0104628 | 0.0545056 | 0.0076278 | 0.6607098 | 36523.62 | 0.2411462 | NA | NA | NA | 0.0216987 | 2513.664 | 0.2094634 | 1.9299725 | 7.892078 | 12.60501 | NA | NA | NA | NA | NA | 0.6233307 | 0.0685092 | 0.0173301 | NA | NA | NA | 0.0252765 | 2.3432735 | 0.6877957 | 0.0371077 | 0.1290979 | 0.1127880 | 0.7102951 |
| 4 | 401 | Winston-Salem | 43.22648 | -0.0193413 | 0.2634850 | 0.0288748 | 0.0478488 | 0.0685591 | 0.4064568 | 30403.07 | 0.4511340 | 9.971601 | 0.3481184 | 0.5376577 | 0.0451547 | 2137.595 | 0.0028042 | -0.0035268 | 7.054153 | 20.89541 | -10.453489 | NA | 0.0418891 | 2523.8348 | 0.1363400 | 0.6357111 | 0.1408433 | 0.8333290 | 0.0239467 | 0.0580380 | 0.0682463 | 0.1903379 | -0.6936674 | 0.3871659 | 0.0499824 | 0.2215046 | 0.1146979 | 0.5695345 |
| 5 | 402 | Martinsville | 53.63568 | 0.0056276 | 0.0069904 | 0.0220931 | 0.0410883 | 0.0172575 | 0.5862997 | 29943.32 | 0.2565523 | 4.503661 | 0.1837683 | 0.7311252 | 0.0266855 | 2756.848 | 0.3987184 | -0.0275933 | 4.544135 | 18.19494 | 5.789057 | -0.0153886 | NA | NA | NA | 0.6121494 | 0.1984681 | 0.5390288 | 0.0377694 | 0.0170314 | 0.0643845 | 0.0236150 | 0.3655432 | 0.8407629 | 0.0321638 | 0.1507503 | 0.1158269 | 0.6865350 |
Nous voyons qu’il y a des données manquantes et que les variables explicatives sont normalement numériques. Nous nous en occuperons plus tard. Pour une explication des variables, considérez le fichier codebook.xls.
select
En général, un jeu de données comporte de nombreuses colonnes. Mais parfois, nous n’avons besoin de considérer que certaines d’entre elles. La logique “select” vous permet de choisir un sous-ensemble des colonnes de données.
Il existe de nombreuses façons de le faire. Ci-dessous la démonstration pour sélectionner une seule colonne “e_rank_b” et une collection de colonnes c(“e_rank_b”, “taxrate”, “crime_violent”). C’est quoi la difference?
# Une seule colonne
facon_1 <- mobility_data$e_rank_b
facon_2 <- mobility_data[,'e_rank_b']
facon_3 <- mobility_data %>% select(e_rank_b)
# Plusieurs
facon_1_mul <- mobility_data[, c("e_rank_b", "taxrate", "crime_violent")]
facon_2_mul <- mobility_data %>% select(e_rank_b, taxrate, crime_violent)
facon_1 %>% head() ## [1] 43.58841 39.64703 46.93892 NA 43.22648 53.63568
facon_1_mul %>% head()## # A tibble: 6 × 3
## e_rank_b taxrate crime_violent
## <dbl> <dbl> <dbl>
## 1 43.6 0.0159 0.0483
## 2 39.6 0.0257 0.0316
## 3 46.9 0.0117 0.0512
## 4 NA 0.0217 0.0371
## 5 43.2 0.0452 0.0500
## 6 53.6 0.0267 0.0322
facon_2 %>% head() ## # A tibble: 6 × 1
## e_rank_b
## <dbl>
## 1 43.6
## 2 39.6
## 3 46.9
## 4 NA
## 5 43.2
## 6 53.6
facon_2_mul %>% head() ## # A tibble: 6 × 3
## e_rank_b taxrate crime_violent
## <dbl> <dbl> <dbl>
## 1 43.6 0.0159 0.0483
## 2 39.6 0.0257 0.0316
## 3 46.9 0.0117 0.0512
## 4 NA 0.0217 0.0371
## 5 43.2 0.0452 0.0500
## 6 53.6 0.0267 0.0322
if-condition (filter)
Les commandes de “select” agissent sur les colonnes, les équivalents pour les lignes sont les “if-conditions” dans R, on appel ca les conditions “filter”. En reprenant l’exemple ci-dessus. Nous ne voulons que les observations où le e_rank_b est supérieur à 40. Quelle différence peut-on constater entre les deux méthodes ?
facon_1 <- mobility_data[mobility_data$e_rank_b > 40,
c("e_rank_b", "taxrate", "crime_violent")]
facon_2 <- mobility_data %>%
select(e_rank_b, taxrate, crime_violent) %>%
filter(e_rank_b > 40)
facon_1 %>% head()## # A tibble: 6 × 3
## e_rank_b taxrate crime_violent
## <dbl> <dbl> <dbl>
## 1 43.6 0.0159 0.0483
## 2 46.9 0.0117 0.0512
## 3 NA NA NA
## 4 43.2 0.0452 0.0500
## 5 53.6 0.0267 0.0322
## 6 45.6 0.0412 0.0267
facon_2 %>% head()## # A tibble: 6 × 3
## e_rank_b taxrate crime_violent
## <dbl> <dbl> <dbl>
## 1 43.6 0.0159 0.0483
## 2 46.9 0.0117 0.0512
## 3 43.2 0.0452 0.0500
## 4 53.6 0.0267 0.0322
## 5 45.6 0.0412 0.0267
## 6 48.8 0.0222 0.0286
groupby
Parfois, nous sommes intéressés par des données agrégées. Avec la commande group_by, nous pouvons facilement créer des statistiques sommaires pour des sous-groupes de données.
Il existe des façons de faire exactement cela dans base R, mais ici nous nous concentrons sur la logique de dplyr, qui est plus intuitive.
# Calculer la moyenne de "e_rank_b" pour chaque ville
group_stat <- mobility_data %>% # Choisir jeu des données
select(czname, e_rank_b) %>% # Choisir les colonnes
group_by(czname) %>% # Specifier les groups
summarise(
mean_per_city = mean(e_rank_b, na.rm=T) # Specifier la statistique d'aggreagtion
)
group_stat %>% head()## # A tibble: 6 × 2
## czname mean_per_city
## <chr> <dbl>
## 1 Aberdeen 45.6
## 2 Aiken 38.5
## 3 Ainsworth 48.6
## 4 Albany 52.9
## 5 Alexandria 46.6
## 6 Allentown 61.3
(mutate)
Souvent, nous devons transformer les variables pour obtenir des informations significatives de notre analyse. Cela peut être fait facilement avec la logique “mutate”. L’utilisation de condtions if-else pour créer de nouvelles variables est extrêmement utile. Par exemple, nous pouvons diviser la mobilité sociale (e_rank_b) en deux catégories : faible et élevée.
data_tmp <- mobility_data
data_tmp <- data_tmp[!is.na(data_tmp$e_rank_b), ]
data_tmp$binary_mobility <- ifelse(data_tmp$e_rank_b > 45, 'High', 'Low')
ex_mutate_2 <- mobility_data %>%
filter(!is.na(e_rank_b)) %>% # PQ?
mutate(binary_mobility = if_else(e_rank_b > 45, 'High', 'Low')) %>%
select(binary_mobility, e_rank_b)
data_tmp[,c('binary_mobility', 'e_rank_b')] %>% head()## # A tibble: 6 × 2
## binary_mobility e_rank_b
## <chr> <dbl>
## 1 Low 43.6
## 2 Low 39.6
## 3 High 46.9
## 4 Low 43.2
## 5 High 53.6
## 6 Low 38.1
ex_mutate_2 %>% head()## # A tibble: 6 × 2
## binary_mobility e_rank_b
## <chr> <dbl>
## 1 Low 43.6
## 2 Low 39.6
## 3 High 46.9
## 4 Low 43.2
## 5 High 53.6
## 6 Low 38.1
2.2 Exercises
- Transformer la variable “e_rank_b” en “inférieure à 50” et “supérieure à 50”. Calculez la moyenne de “crime_violent” pour chaque groupe. Quelles sont vos conclusions?
- Transformer la variable “frac_middleclass” en “inférieure à 0.5” et “supérieure à 0.5”. Pareil pour la variable “eitc_exposure”. Calculer la moyenne du taux d’abandon scolaire au lycée “dropout_r” pour chaque combinaison des deux nouvelles variables. Conseil : utilisez group_by
- Proposer une façon de traiter les valeurs “NA”. Pourquoi pensez-vous que c’est une bonne idée ?
3 Analyse Visuelle
Il est toujours bon d’avoir une compréhension intuitive des données en utilisant des graphiques. Ici, nous essaierons de réaliser la plupart des graphiques à la main (il existe des logiciels qui peuvent faire le travail pour vous, mais il est préférable de le comprendre d’abord complètement).
3.1 Plots
Les plots/figures peuvent facilement être intégrées dans markdown. Tout comme une sortie de tableau, on peut les structurer légèrement. Par exemple, nous pouvons ajouter une légende ou ajouter plusieurs graphiques côte à côte.
Par exemple, commençons par un histogramme:
# Histogram simple
hist(mobility_data$e_rank_b)?hist
?histplot(density(data_tmp$e_rank_b)) # Pourquoi changer les données?
plot(density(data_tmp$e_rank_b), col='red')Il existe aussi un package sympa de l’écosystème tidyverse (ggplot2)
plot_1 <- mobility_data %>%
ggplot() +
geom_histogram(aes(x=e_rank_b))
plot_2 <- mobility_data %>%
ggplot() +
geom_histogram(aes(x=crime_violent))
grid.arrange(plot_1, plot_2, ncol=2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 19 rows containing non-finite values (stat_bin).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 17 rows containing non-finite values (stat_bin).
Nous pouvons également construire des figures avec plusieurs informations. Par exemple, un histogramme avec une ligne de densité. (Attention toutefois aux axes!)
# Base R
hist(data_tmp$e_rank_b, probability = T)
lines(density(data_tmp$e_rank_b), col='red')# GGplot
data_tmp %>%
ggplot() +
geom_histogram(aes(x=e_rank_b, y=..density..)) +
geom_density(aes(x=e_rank_b), color='red')## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
3.2 Exercises
Conseil : si vous n’arrivez pas à trouver la bonne commande pour une figure, les solutions, utilisez google!
Montrer la relation entre le taux d’obtention d’un diplôme universitaire et la mobilité socialeà l’aide d’un diagramme de dispersion (scatterplot) et d’une ligne de régression. Conseil: vous pouvez trouver une solution si vous googlez “r plot with regression line”
Existe-t-il une relation entre les dépenses des gouvernements locaux par habitant (
subcty_total_expenditure_pc) et le taux d’imposition (taxrate)? Pouvez-vous utiliser les données directement ? Conseil : pour mettre les données à l’échelle, vous pouvez toujours essayer d’utiliser log, sqrt, etc.Montrez graphiquement qu’il existe une ville où la criminalité violente est exceptionnellement élevée (‘crime_violent’). Que pouvez-vous faire avec de telles données (
outlierou pas)?
4 Analyse descriptive
L’analyse graphique est utile mais se complique avec plus que deux variables. En général, nous évaluons également les hypothèses que nous avons dérivées des graphiques à l’aide d’une sorte de test.
4.1 Tables
Les tableaux peuvent facilement être incorporés dans R. Vous pouvez soit les imprimer directement, soit utiliser kable (du paquet knitr - plus jolie). Bien sûr, nous pouvons également combiner ce que nous avons vu dans les sections précédentes (en particulier la manipulation des données)
# Impression de base
mobility_data %>% head()## # A tibble: 6 × 39
## index cz czname e_rank_b cs_race_bla cs_race_theil_20… cs00_seg_inc
## <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 0 100 Johnson City 43.6 0.0694 0.0577 0.0409
## 2 1 200 Morristown 39.6 0.445 0.0750 -0.0259
## 3 2 301 Middlesborough 46.9 0.0151 0.110 0.0954
## 4 3 302 Knoxville NA -0.000411 -0.00446 0.0105
## 5 4 401 Winston-Salem 43.2 -0.0193 0.263 0.0289
## 6 5 402 Martinsville 53.6 0.00563 0.00699 0.0221
## # … with 32 more variables: cs00_seg_inc_pov25 <dbl>, cs00_seg_inc_aff75 <dbl>,
## # frac_traveltime_lt15 <dbl>, hhinc00 <dbl>, gini <dbl>,
## # inc_share_1perc <dbl>, gini99 <dbl>, frac_middleclass <dbl>, taxrate <dbl>,
## # subcty_total_expenditure_pc <dbl>, tax_st_diff_top20 <dbl>,
## # eitc_exposure <dbl>, ccd_exp_tot <dbl>, ccd_pup_tch_ratio <dbl>,
## # score_r <dbl>, dropout_r <dbl>, num_inst_pc <dbl>, tuition <dbl>,
## # gradrate_r <dbl>, cs_labforce <dbl>, cs_elf_ind_man <dbl>, …
# Avec Kable
mobility_data %>% head() %>% kable() # c'est comme kable(head(mobilitiy_data))| index | cz | czname | e_rank_b | cs_race_bla | cs_race_theil_2000 | cs00_seg_inc | cs00_seg_inc_pov25 | cs00_seg_inc_aff75 | frac_traveltime_lt15 | hhinc00 | gini | inc_share_1perc | gini99 | frac_middleclass | taxrate | subcty_total_expenditure_pc | tax_st_diff_top20 | eitc_exposure | ccd_exp_tot | ccd_pup_tch_ratio | score_r | dropout_r | num_inst_pc | tuition | gradrate_r | cs_labforce | cs_elf_ind_man | d_tradeusch_pw_1990 | frac_worked1416 | mig_inflow | mig_outflow | cs_born_foreign | scap_ski90pcm | rel_tot | crime_violent | cs_fam_wkidsinglemom | cs_divorced | cs_married |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 100 | Johnson City | 43.58841 | 0.0693926 | 0.0576880 | 0.0408645 | -0.0171138 | 0.0266398 | 0.3844612 | 38209.54 | 0.4193436 | 5.695437 | 0.2528637 | 0.5840611 | 0.0158913 | 1404.822 | 1.0794937 | 0.0000780 | 6.577284 | 17.96812 | NA | NA | 0.0554321 | 1032.0113 | -0.2210098 | 0.6529477 | 0.0600338 | NA | 0.0148774 | 0.0417411 | 0.0559070 | 0.1604878 | NA | 0.5358987 | 0.0482720 | 0.2421846 | 0.1392865 | 0.5491035 |
| 1 | 200 | Morristown | 39.64703 | 0.4453178 | 0.0749721 | -0.0258754 | 0.0352996 | 0.0089049 | 0.3543061 | 25123.34 | 0.4867183 | 11.816536 | 0.3926104 | 0.3559214 | 0.0257373 | 1568.710 | 0.0198662 | -0.0128659 | 4.405095 | 17.79943 | -16.015344 | 0.0561887 | 0.0379247 | 970.0455 | -0.0416078 | 0.5273348 | 0.2086723 | 0.7670213 | 0.0223866 | 0.0181590 | 0.0137850 | 0.0232165 | -0.6127989 | 0.6469429 | 0.0316044 | 0.3499764 | 0.1039057 | 0.5426201 |
| 2 | 301 | Middlesborough | 46.93892 | 0.0151223 | 0.1100580 | 0.0954185 | 0.0139200 | 0.0527287 | 0.4793166 | 32701.43 | 0.3071970 | 11.810008 | 0.2077348 | 0.6675747 | 0.0116976 | 3115.946 | -0.0046168 | 0.0343010 | 6.609784 | 17.85102 | -9.137633 | 0.0246227 | 0.0634999 | 1395.0190 | -0.1716805 | 0.7066224 | 0.0362908 | 0.0274560 | 0.0274638 | 0.0309692 | 0.0421808 | 0.1318125 | -2.0635338 | 0.4222386 | 0.0511852 | 0.1831544 | 0.1402021 | 0.6224286 |
| 3 | 302 | Knoxville | NA | -0.0004112 | -0.0044561 | 0.0104628 | 0.0545056 | 0.0076278 | 0.6607098 | 36523.62 | 0.2411462 | NA | NA | NA | 0.0216987 | 2513.664 | 0.2094634 | 1.9299725 | 7.892078 | 12.60501 | NA | NA | NA | NA | NA | 0.6233307 | 0.0685092 | 0.0173301 | NA | NA | NA | 0.0252765 | 2.3432735 | 0.6877957 | 0.0371077 | 0.1290979 | 0.1127880 | 0.7102951 |
| 4 | 401 | Winston-Salem | 43.22648 | -0.0193413 | 0.2634850 | 0.0288748 | 0.0478488 | 0.0685591 | 0.4064568 | 30403.07 | 0.4511340 | 9.971601 | 0.3481184 | 0.5376577 | 0.0451547 | 2137.595 | 0.0028042 | -0.0035268 | 7.054153 | 20.89541 | -10.453489 | NA | 0.0418891 | 2523.8348 | 0.1363400 | 0.6357111 | 0.1408433 | 0.8333290 | 0.0239467 | 0.0580380 | 0.0682463 | 0.1903379 | -0.6936674 | 0.3871659 | 0.0499824 | 0.2215046 | 0.1146979 | 0.5695345 |
| 5 | 402 | Martinsville | 53.63568 | 0.0056276 | 0.0069904 | 0.0220931 | 0.0410883 | 0.0172575 | 0.5862997 | 29943.32 | 0.2565523 | 4.503661 | 0.1837683 | 0.7311252 | 0.0266855 | 2756.848 | 0.3987184 | -0.0275933 | 4.544135 | 18.19494 | 5.789057 | -0.0153886 | NA | NA | NA | 0.6121494 | 0.1984681 | 0.5390288 | 0.0377694 | 0.0170314 | 0.0643845 | 0.0236150 | 0.3655432 | 0.8407629 | 0.0321638 | 0.1507503 | 0.1158269 | 0.6865350 |
ici nous pouvons manipuler légèrement nos données
# Example avec manipulation des données
# Décrire ce qui se passe ici
mobility_data %>%
filter(!is.na(e_rank_b)) %>%
mutate(decile = ntile(e_rank_b, 10)) %>%
select(decile, e_rank_b, mig_outflow) %>%
group_by(decile) %>%
summarise(average_outflow = mean(mig_outflow)) %>%
kable()| decile | average_outflow |
|---|---|
| 1 | 0.0463487 |
| 2 | 0.0459592 |
| 3 | 0.0457738 |
| 4 | 0.0436270 |
| 5 | 0.0413216 |
| 6 | 0.0447277 |
| 7 | 0.0432519 |
| 8 | 0.0408539 |
| 9 | 0.0399971 |
| 10 | NA |
4.2 Distributions en R
L’un des exercices les plus importants lors de la modélisation statistique consiste à vérifier les hypothèses du modèle. La logique est toujours la même.
q-distribution (quantiles, ex: qnorm, qpois, qgamma etc.) Quel est le 99ème quantile d’une gamma avec shape=2 et rate=1?
qgamma(0.99, shape=2,rate=1)## [1] 6.638352
d-distribution (densité, ex: dnorm, dpois) Tracez visuellement le 99ème quantile d’une gamma avec shape=2 et rate=1?
gamma_dist <- dgamma(seq(0,8,0.1), shape=2,rate=1)
plot(x=seq(0,8,0.1),
y=gamma_dist,
type='l')
abline(v=qgamma(0.99, shape=2,rate=1),
col='red', lty=2)p-distribution (densité cumulative, ex: pnorm, ppois)
Supposons que \(X \sim Gamma(2,1)\). Quel pourcentage de la masse totale se trouve à droite de 5 ?
1 - pgamma(5, shape=2, rate=1)## [1] 0.04042768
r-distribution (échantillon, ex: rnorm, rpois)
Simulez 100 observations à partir d’un gamma avec shape=2, rate=1 et comparez-les à la distribution théorique.
set.seed(42)
sims_gamma <- rgamma(100, shape=2, rate=1)
plot(density(sims_gamma), ylim=c(0,0.38),
xlim=c(-1, 7),
main='Comparaison empirique vs. theoretique')
lines(x=seq(0,8,0.1),
y=gamma_dist, col='red')4.2.1 Q-Q Plot
Pour comparer un ajustement avec une distribution théorique, on peut utiliser un diagramme quantile-quantile (qq-plot). Par exemple, on suppose souvent que le log revenue est distribué selon une normale. Comparons l’ajustement théorique et l’ajustement réel.
# A la main
# Sans transformation
inc <- (as.numeric(mobility_data$hhinc00))
inc <- (inc - mean(inc)) / sd(inc)
quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(inc, seq(0,1,0.001))
plot(quantiles_norm,
quantiles_empirical,
main='Normal Q-Q Plot',
ylab='Sample Quantiles',
xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")# Avec transformation
# A la main
# Sans transformation
log_inc <- log(as.numeric(mobility_data$hhinc00))
log_inc <- (log_inc - mean(log_inc)) / sd(log_inc)
quantiles_norm <- qnorm(seq(0,1,0.001))
quantiles_empirical <- quantile(log_inc, seq(0,1,0.001))
plot(quantiles_norm,
quantiles_empirical,
main='Normal Q-Q Plot (Log transformed)',
ylab='Sample Quantiles',
xlab='Theoretical Quantiles')
abline(a = 0, b = 1, col = "red")qqnorm(log(as.numeric(mobility_data$hhinc00)))
qqline(log(as.numeric(mobility_data$hhinc00)),
col = 2,lwd=2,lty=2)4.3 Exercises
Supposons que \(X \sim \mathcal{N}(0,10)\). Echantillonnez 10 observations et calculez la moyenne, que nous appelons Y. Répétez l’opération ci-dessus 1000 fois et calculez la moyenne et l’écart-type de Y. A quoi vous attendez-vous et pourquoi, qu’observez-vous ?
Créez un tableau ou qui montre le nombre d’observations manquantes pour chaque variable (colonne dans l’ensemble de données).
Analysez les variables
tuitionetnum_inst_pcet discutez de ce que vous voyez. Quelles sont vos conclusions ? S’agit-il d’un problème?Créez un tableau montrant le coefficient de corrélation entre la mobilité sociale
e_rank_bet toutes les autres variables. Évaluez vos résultats.Séparez la variable en “faible”, “moyen” et “élevé”. Justifiez votre choix. Comparez chaque niveau avec les variables
ginietcs_born_foreign. Interprétez et créez un tableau.
5 Example guidée
Dans cet exemple guidé, nous allons à nouveau travailler avec les données de Raj Chetty, Nathaniel Hendren, Patrick Kline et Emmanuel Saez.
Notre objectif est de mieux comprendre la mobilité intragénérationnelle. On dit que la mobilité intergénérationnelle est élevée si le revenu des parents (au cours de leur vie) n’est pas très lié au revenu de leurs enfants. Pour les exercices, vous pouvez utiliser n’importe quel paquet que vous jugez nécessaire mais vous devez justifier vos résultats.
Pour cela, résolvez les exercices suivants :
Chargez les données, regardez les types de données, est-ce que cela a du sens ?
La variable d’intérêt est
e_rank_b, renommez cette colonne. Que se passe-t-il si vous utilisez un espace ” ” dans le nom ? Est-ce une bonne idée ?Visualiser la distribution de la mobilité intergénérationnelle. Donnez à vos graphiques un titre approprié et des étiquettes d’axe correctes.
Examinez la relation entre la mobilité intergénérationnelle et la taille du secteur manufacturier. Quelle est votre conclusion ?
Découpez la variable de la taille du secteur manufacturier en trois composantes (low-mid-high). Calculez les moyennes par groupe pour trois autres variables que vous jugez intéressantes.
Ajustez une régression qui décrit l’équation suivante \(\mathbb{E}[y| \text{manufacturing}]=\hat{\beta}_0 + \hat{\beta}_1*\text{manuf_1} + \hat{\beta}_2*\text{manuf_2}\). Pourquoi exclure la dernière composante?
L’existence d’une classe moyenne importante est souvent considérée comme une bonne chose. Voir par exemple cet article, analysons cela plus en détail.
Réalisez un nuage de points (“scatterplot”) avec la mobilité intergénérationnelle sur l’axe des y et la fraction de la classe moyenne sur l’axe des x. Quelles sont vos conclusions ?
Exécutez un modèle linéaire qui explique la mobilité intergénérationnelle à l’aide de la fraction de la classe moyenne.
Nous voulons rendre le modèle légèrement plus flexible. Effectuez 10 régressions différentes qui correspondent à des polynômes de plus en plus complexes de la fraction de la classe moyenne.
Défi : Visualisez les différents résultats. Qu’arrive-t-il à la valeur \(R^2\)? Est-ce une bonne chose ?
5.1 Charger et preparer les données
# Charger l'essentiel
library(tidyverse)
library(knitr)
library(gridExtra)
library(haven)
library(readxl)
mobility_data <- read_dta("data/mobility_across_cz.dta")
codebook <- read_excel("data/codebook.xls")
mobility_data %>%
head() %>%
kable()| index | cz | czname | e_rank_b | cs_race_bla | cs_race_theil_2000 | cs00_seg_inc | cs00_seg_inc_pov25 | cs00_seg_inc_aff75 | frac_traveltime_lt15 | hhinc00 | gini | inc_share_1perc | gini99 | frac_middleclass | taxrate | subcty_total_expenditure_pc | tax_st_diff_top20 | eitc_exposure | ccd_exp_tot | ccd_pup_tch_ratio | score_r | dropout_r | num_inst_pc | tuition | gradrate_r | cs_labforce | cs_elf_ind_man | d_tradeusch_pw_1990 | frac_worked1416 | mig_inflow | mig_outflow | cs_born_foreign | scap_ski90pcm | rel_tot | crime_violent | cs_fam_wkidsinglemom | cs_divorced | cs_married |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 100 | Johnson City | 43.58841 | 0.0693926 | 0.0576880 | 0.0408645 | -0.0171138 | 0.0266398 | 0.3844612 | 38209.54 | 0.4193436 | 5.695437 | 0.2528637 | 0.5840611 | 0.0158913 | 1404.822 | 1.0794937 | 0.0000780 | 6.577284 | 17.96812 | NA | NA | 0.0554321 | 1032.0113 | -0.2210098 | 0.6529477 | 0.0600338 | NA | 0.0148774 | 0.0417411 | 0.0559070 | 0.1604878 | NA | 0.5358987 | 0.0482720 | 0.2421846 | 0.1392865 | 0.5491035 |
| 1 | 200 | Morristown | 39.64703 | 0.4453178 | 0.0749721 | -0.0258754 | 0.0352996 | 0.0089049 | 0.3543061 | 25123.34 | 0.4867183 | 11.816536 | 0.3926104 | 0.3559214 | 0.0257373 | 1568.710 | 0.0198662 | -0.0128659 | 4.405095 | 17.79943 | -16.015344 | 0.0561887 | 0.0379247 | 970.0455 | -0.0416078 | 0.5273348 | 0.2086723 | 0.7670213 | 0.0223866 | 0.0181590 | 0.0137850 | 0.0232165 | -0.6127989 | 0.6469429 | 0.0316044 | 0.3499764 | 0.1039057 | 0.5426201 |
| 2 | 301 | Middlesborough | 46.93892 | 0.0151223 | 0.1100580 | 0.0954185 | 0.0139200 | 0.0527287 | 0.4793166 | 32701.43 | 0.3071970 | 11.810008 | 0.2077348 | 0.6675747 | 0.0116976 | 3115.946 | -0.0046168 | 0.0343010 | 6.609784 | 17.85102 | -9.137633 | 0.0246227 | 0.0634999 | 1395.0190 | -0.1716805 | 0.7066224 | 0.0362908 | 0.0274560 | 0.0274638 | 0.0309692 | 0.0421808 | 0.1318125 | -2.0635338 | 0.4222386 | 0.0511852 | 0.1831544 | 0.1402021 | 0.6224286 |
| 3 | 302 | Knoxville | NA | -0.0004112 | -0.0044561 | 0.0104628 | 0.0545056 | 0.0076278 | 0.6607098 | 36523.62 | 0.2411462 | NA | NA | NA | 0.0216987 | 2513.664 | 0.2094634 | 1.9299725 | 7.892078 | 12.60501 | NA | NA | NA | NA | NA | 0.6233307 | 0.0685092 | 0.0173301 | NA | NA | NA | 0.0252765 | 2.3432735 | 0.6877957 | 0.0371077 | 0.1290979 | 0.1127880 | 0.7102951 |
| 4 | 401 | Winston-Salem | 43.22648 | -0.0193413 | 0.2634850 | 0.0288748 | 0.0478488 | 0.0685591 | 0.4064568 | 30403.07 | 0.4511340 | 9.971601 | 0.3481184 | 0.5376577 | 0.0451547 | 2137.595 | 0.0028042 | -0.0035268 | 7.054153 | 20.89541 | -10.453489 | NA | 0.0418891 | 2523.8348 | 0.1363400 | 0.6357111 | 0.1408433 | 0.8333290 | 0.0239467 | 0.0580380 | 0.0682463 | 0.1903379 | -0.6936674 | 0.3871659 | 0.0499824 | 0.2215046 | 0.1146979 | 0.5695345 |
| 5 | 402 | Martinsville | 53.63568 | 0.0056276 | 0.0069904 | 0.0220931 | 0.0410883 | 0.0172575 | 0.5862997 | 29943.32 | 0.2565523 | 4.503661 | 0.1837683 | 0.7311252 | 0.0266855 | 2756.848 | 0.3987184 | -0.0275933 | 4.544135 | 18.19494 | 5.789057 | -0.0153886 | NA | NA | NA | 0.6121494 | 0.1984681 | 0.5390288 | 0.0377694 | 0.0170314 | 0.0643845 | 0.0236150 | 0.3655432 | 0.8407629 | 0.0321638 | 0.1507503 | 0.1158269 | 0.6865350 |
str(mobility_data)## tibble [495 × 39] (S3: tbl_df/tbl/data.frame)
## $ index : num [1:495] 0 1 2 3 4 5 6 7 8 9 ...
## ..- attr(*, "format.stata")= chr "%12.0g"
## $ cz : num [1:495] 100 200 301 302 401 402 500 601 602 700 ...
## ..- attr(*, "format.stata")= chr "%12.0g"
## $ czname : chr [1:495] "Johnson City" "Morristown" "Middlesborough" "Knoxville" ...
## ..- attr(*, "format.stata")= chr "%22s"
## $ e_rank_b : num [1:495] 43.6 39.6 46.9 NA 43.2 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_race_bla : num [1:495] 0.069393 0.445318 0.015122 -0.000411 -0.019341 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_race_theil_2000 : num [1:495] 0.05769 0.07497 0.11006 -0.00446 0.26349 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs00_seg_inc : num [1:495] 0.0409 -0.0259 0.0954 0.0105 0.0289 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs00_seg_inc_pov25 : num [1:495] -0.0171 0.0353 0.0139 0.0545 0.0478 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs00_seg_inc_aff75 : num [1:495] 0.02664 0.0089 0.05273 0.00763 0.06856 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ frac_traveltime_lt15 : num [1:495] 0.384 0.354 0.479 0.661 0.406 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ hhinc00 : num [1:495] 38210 25123 32701 36524 30403 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ gini : num [1:495] 0.419 0.487 0.307 0.241 0.451 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ inc_share_1perc : num [1:495] 5.7 11.82 11.81 NA 9.97 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ gini99 : num [1:495] 0.253 0.393 0.208 NA 0.348 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ frac_middleclass : num [1:495] 0.584 0.356 0.668 NA 0.538 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ taxrate : num [1:495] 0.0159 0.0257 0.0117 0.0217 0.0452 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ subcty_total_expenditure_pc: num [1:495] 1405 1569 3116 2514 2138 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ tax_st_diff_top20 : num [1:495] 1.07949 0.01987 -0.00462 0.20946 0.0028 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ eitc_exposure : num [1:495] 0.000078 -0.012866 0.034301 1.929972 -0.003527 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ ccd_exp_tot : num [1:495] 6.58 4.41 6.61 7.89 7.05 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ ccd_pup_tch_ratio : num [1:495] 18 17.8 17.9 12.6 20.9 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ score_r : num [1:495] NA -16.02 -9.14 NA -10.45 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ dropout_r : num [1:495] NA 0.0562 0.0246 NA NA ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ num_inst_pc : num [1:495] 0.0554 0.0379 0.0635 NA 0.0419 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ tuition : num [1:495] 1032 970 1395 NA 2524 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ gradrate_r : num [1:495] -0.221 -0.0416 -0.1717 NA 0.1363 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_labforce : num [1:495] 0.653 0.527 0.707 0.623 0.636 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_elf_ind_man : num [1:495] 0.06 0.2087 0.0363 0.0685 0.1408 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ d_tradeusch_pw_1990 : num [1:495] NA 0.767 0.0275 0.0173 0.8333 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ frac_worked1416 : num [1:495] 0.0149 0.0224 0.0275 NA 0.0239 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ mig_inflow : num [1:495] 0.0417 0.0182 0.031 NA 0.058 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ mig_outflow : num [1:495] 0.0559 0.0138 0.0422 NA 0.0682 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_born_foreign : num [1:495] 0.1605 0.0232 0.1318 0.0253 0.1903 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ scap_ski90pcm : num [1:495] NA -0.613 -2.064 2.343 -0.694 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ rel_tot : num [1:495] 0.536 0.647 0.422 0.688 0.387 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ crime_violent : num [1:495] 0.0483 0.0316 0.0512 0.0371 0.05 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_fam_wkidsinglemom : num [1:495] 0.242 0.35 0.183 0.129 0.222 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_divorced : num [1:495] 0.139 0.104 0.14 0.113 0.115 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ cs_married : num [1:495] 0.549 0.543 0.622 0.71 0.57 ...
## ..- attr(*, "format.stata")= chr "%10.0g"
Le format est celui de stata, mais la plupart des variables sont numériques, cela semble correct.
# Verifier données manquantes
for (column_ in names(mobility_data)){
print(
paste(
column_,
sum(is.na(mobility_data[, column_])),
sep=': '
)
)
}## [1] "index: 0"
## [1] "cz: 0"
## [1] "czname: 0"
## [1] "e_rank_b: 19"
## [1] "cs_race_bla: 0"
## [1] "cs_race_theil_2000: 0"
## [1] "cs00_seg_inc: 0"
## [1] "cs00_seg_inc_pov25: 0"
## [1] "cs00_seg_inc_aff75: 0"
## [1] "frac_traveltime_lt15: 0"
## [1] "hhinc00: 0"
## [1] "gini: 0"
## [1] "inc_share_1perc: 19"
## [1] "gini99: 19"
## [1] "frac_middleclass: 19"
## [1] "taxrate: 1"
## [1] "subcty_total_expenditure_pc: 2"
## [1] "tax_st_diff_top20: 0"
## [1] "eitc_exposure: 0"
## [1] "ccd_exp_tot: 7"
## [1] "ccd_pup_tch_ratio: 20"
## [1] "score_r: 21"
## [1] "dropout_r: 89"
## [1] "num_inst_pc: 100"
## [1] "tuition: 104"
## [1] "gradrate_r: 107"
## [1] "cs_labforce: 0"
## [1] "cs_elf_ind_man: 0"
## [1] "d_tradeusch_pw_1990: 14"
## [1] "frac_worked1416: 19"
## [1] "mig_inflow: 10"
## [1] "mig_outflow: 10"
## [1] "cs_born_foreign: 0"
## [1] "scap_ski90pcm: 14"
## [1] "rel_tot: 0"
## [1] "crime_violent: 17"
## [1] "cs_fam_wkidsinglemom: 0"
## [1] "cs_divorced: 0"
## [1] "cs_married: 0"
Le nombre de points de données manquants ne semble pas non plus trop inquiétant. Il pourrait y avoir une certaine corrélation entre les points manquants mais nous pouvons ignorer cela pour le moment.
5.2 Renommer une variable
# renommer avec espace
print(names(mobility_data))## [1] "index" "cz"
## [3] "czname" "e_rank_b"
## [5] "cs_race_bla" "cs_race_theil_2000"
## [7] "cs00_seg_inc" "cs00_seg_inc_pov25"
## [9] "cs00_seg_inc_aff75" "frac_traveltime_lt15"
## [11] "hhinc00" "gini"
## [13] "inc_share_1perc" "gini99"
## [15] "frac_middleclass" "taxrate"
## [17] "subcty_total_expenditure_pc" "tax_st_diff_top20"
## [19] "eitc_exposure" "ccd_exp_tot"
## [21] "ccd_pup_tch_ratio" "score_r"
## [23] "dropout_r" "num_inst_pc"
## [25] "tuition" "gradrate_r"
## [27] "cs_labforce" "cs_elf_ind_man"
## [29] "d_tradeusch_pw_1990" "frac_worked1416"
## [31] "mig_inflow" "mig_outflow"
## [33] "cs_born_foreign" "scap_ski90pcm"
## [35] "rel_tot" "crime_violent"
## [37] "cs_fam_wkidsinglemom" "cs_divorced"
## [39] "cs_married"
names(mobility_data)[4] <- "mobilite intergen"
# mobility_data$`mobilite intergen` -> On a toujours besoin des symboles `. Cela semble encombrant.
# Essayez d'utiliser des tirets bas (_) ou des points (.)
names(mobility_data)[4] <- "mobilite_intergen"
# mobility_data$mobilite_intergen -> C'est mieux5.3 Analyse visuelle de la mobilité
Nous allons analyser la densité en utilisant le lissage par noyau, un histogramme et la fonction de distribution cumulative.
hist(mobility_data$mobilite_intergen,
main='Intergenerational Mobility Scores',
xlab='Value',
ylab='Density',
probability = T,
breaks = seq(27, 67, 3))
lines(density(mobility_data$mobilite_intergen[!is.na(mobility_data$mobilite_intergen)]),
col='red',
lwd=2,
lty=2)
legend(55,0.07,
legend=c("Density"),
col=c("red"),
lty=2,
cex=0.8)quantiles_ <- quantile(mobility_data$mobilite_intergen,
probs = seq(0,1,0.01),
na.rm = T)
plot(x=quantiles_,
y=seq(0,1,0.01),
type='l',
lwd=2,
main='ECDF',
ylab='',
xlab='Value for Mobility')5.4 Analyse visuelle: Manufacturing vs Mobilité
plot(
mobility_data$cs_elf_ind_man,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
) # Model simple pour la relation lineaire
summary(lm(mobilite_intergen ~ cs_elf_ind_man, data=mobility_data))##
## Call:
## lm(formula = mobilite_intergen ~ cs_elf_ind_man, data = mobility_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5226 -3.6620 -0.8982 3.1683 17.9719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.9024 0.5867 79.946 < 2e-16 ***
## cs_elf_ind_man -17.7064 3.1306 -5.656 2.68e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.482 on 474 degrees of freedom
## (19 observations deleted due to missingness)
## Multiple R-squared: 0.06322, Adjusted R-squared: 0.06125
## F-statistic: 31.99 on 1 and 474 DF, p-value: 2.683e-08
plot(
mobility_data$cs_elf_ind_man,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
)
abline(a=46.9024,
b=-17.7064,
col='red',
lwd=2,
lty=2)
legend(0.35,60,
legend=c("Linear Fit"),
col=c("red"),
lty=2,
cex=0.8)5.5 Analyse descriptif (secteur manufacturier)
Nous pouvons trouver de bons seuils en regardant le graphique ci-dessus ou en utilisant les quantiles. Voici l’approche “visuelle”. Lorsque nous coupons une variable, nous pouvons introduire des non-linéarités dans la relation (voir l’exercice suivant)
mobility_data$manuf_cut <- cut(mobility_data$cs_elf_ind_man,
breaks = c(0,0.1,0.3, 1),
labels = c('low', 'mid', 'high'))
mobility_data %>%
group_by(manuf_cut) %>%
summarise(
mean_labforce = mean(cs_labforce, na.rm=T),
mean_growth_imports = mean(d_tradeusch_pw_1990, na.rm=T),
mean_mig_inflow = mean(mig_outflow, na.rm=T)
)## # A tibble: 3 × 4
## manuf_cut mean_labforce mean_growth_imports mean_mig_inflow
## <fct> <dbl> <dbl> <dbl>
## 1 low 0.633 0.289 0.0458
## 2 mid 0.641 1.30 0.0427
## 3 high 0.667 2.73 0.0399
5.6 Model simple
model_lineaire_simple = lm(mobilite_intergen ~ manuf_cut,
data=mobility_data)
summary(model_lineaire_simple)##
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut, data = mobility_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.271 -3.847 -0.658 2.952 17.605
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.9575 0.5388 85.293 < 2e-16 ***
## manuf_cutmid -2.3677 0.6141 -3.856 0.000131 ***
## manuf_cuthigh -6.3979 1.2444 -5.142 4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.495 on 473 degrees of freedom
## (19 observations deleted due to missingness)
## Multiple R-squared: 0.06088, Adjusted R-squared: 0.05691
## F-statistic: 15.33 on 2 and 473 DF, p-value: 3.534e-07
plot(
mobility_data$cs_elf_ind_man,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
)
x <- seq(0, 0.45, 0.01)
fx <- 45.9575 + (x > 0 & x <=0.1) * 0 +
(x > 0.1 & x <= 0.3) * -2.3677 +
(x > 0.3 & x <= 2.119) * -6.3979
lines(x, fx,
col='red',
lty=2,
lwd=2)plot(
mobility_data$cs_elf_ind_man,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
)
predictions_simple <- predict(model_lineaire_simple,
mobility_data)
points(mobility_data$cs_elf_ind_man,
predictions_simple,
col='red')model_lineaire_interaction = lm(mobilite_intergen ~ manuf_cut + cs_elf_ind_man + manuf_cut*cs_elf_ind_man,
data=mobility_data)
summary(model_lineaire_interaction)##
## Call:
## lm(formula = mobilite_intergen ~ manuf_cut + cs_elf_ind_man +
## manuf_cut * cs_elf_ind_man, data = mobility_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.9773 -3.6249 -0.6577 2.7890 17.3450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.921 1.983 21.640 <2e-16 ***
## manuf_cutmid 3.175 2.252 1.410 0.1592
## manuf_cuthigh -5.908 10.728 -0.551 0.5821
## cs_elf_ind_man 44.485 27.973 1.590 0.1124
## manuf_cutmid:cs_elf_ind_man -57.912 28.507 -2.031 0.0428 *
## manuf_cuthigh:cs_elf_ind_man -37.336 40.598 -0.920 0.3582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.463 on 470 degrees of freedom
## (19 observations deleted due to missingness)
## Multiple R-squared: 0.07768, Adjusted R-squared: 0.06787
## F-statistic: 7.917 on 5 and 470 DF, p-value: 3.562e-07
plot(
mobility_data$cs_elf_ind_man,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
)
predictions_interact <- predict(model_lineaire_interaction,
mobility_data)
points(mobility_data$cs_elf_ind_man,
predictions_interact,
col='red')5.7 Model “Classe moyenne”
plot(x=mobility_data$frac_middleclass,
y=mobility_data$mobilite_intergen,
main='Middle Class and Intergen. Mobility',
ylab='Intergen. Mobility',
xlab='Fraction Middle class')model_middle_class <- lm(mobilite_intergen ~ frac_middleclass,
data=mobility_data)
summary(model_middle_class)##
## Call:
## lm(formula = mobilite_intergen ~ frac_middleclass, data = mobility_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9940 -2.9766 -0.5619 2.7327 15.9744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.986 1.381 13.74 <2e-16 ***
## frac_middleclass 45.411 2.491 18.23 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.343 on 474 degrees of freedom
## (19 observations deleted due to missingness)
## Multiple R-squared: 0.4121, Adjusted R-squared: 0.4109
## F-statistic: 332.3 on 1 and 474 DF, p-value: < 2.2e-16
plot(
mobility_data$frac_middleclass,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
)
predictions_middleclass <- predict(model_middle_class,
mobility_data)
points(mobility_data$frac_middleclass,
predictions_middleclass,
col='red')# Analyse des residues
residuals_ <- mobility_data$mobilite_intergen - predictions_middleclass
plot(mobility_data$frac_middleclass,
residuals_,
ylab='Residuals',
xlab='Middle class')mobility_data_nomissings <- mobility_data[!is.na(mobility_data$frac_middleclass), ]
model_middle_class_poly <- lm(mobilite_intergen ~ poly(frac_middleclass, 2),
data=mobility_data_nomissings)
summary(model_middle_class_poly)##
## Call:
## lm(formula = mobilite_intergen ~ poly(frac_middleclass, 2), data = mobility_data_nomissings)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0057 -2.9207 -0.3494 2.4430 15.5096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.9039 0.1881 233.464 < 2e-16 ***
## poly(frac_middleclass, 2)1 79.1684 4.1029 19.296 < 2e-16 ***
## poly(frac_middleclass, 2)2 31.2715 4.1029 7.622 1.38e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.103 on 473 degrees of freedom
## Multiple R-squared: 0.4764, Adjusted R-squared: 0.4742
## F-statistic: 215.2 on 2 and 473 DF, p-value: < 2.2e-16
plot(
mobility_data$frac_middleclass,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
)
predictions_middleclass_poly <- predict(model_middle_class_poly,
mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])
lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
predictions_middleclass_poly,
col='red')5.8 Challenge
# List pour les resultats
results_list <- list()
# Loop
for (length_ in seq(1,10,1)){
reg_ <- lm(mobilite_intergen ~ poly(frac_middleclass, length_),
data=mobility_data_nomissings)
print(paste('R-Squared: ', summary(reg_)$r.squared))
predictions_ <- predict(reg_,
mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ])
results_list[[length_]] <- predictions_
}## [1] "R-Squared: 0.412134418340672"
## [1] "R-Squared: 0.476437415437812"
## [1] "R-Squared: 0.477700187650047"
## [1] "R-Squared: 0.477727823275877"
## [1] "R-Squared: 0.478923076244006"
## [1] "R-Squared: 0.48017980394324"
## [1] "R-Squared: 0.48426795442976"
## [1] "R-Squared: 0.487546350228303"
## [1] "R-Squared: 0.487606950061351"
## [1] "R-Squared: 0.489193655479349"
plot(
mobility_data$frac_middleclass,
mobility_data$mobilite_intergen,
main='% Manufacturing vs Mobility',
ylab='Intergen. Mob',
xlab='Fraction Working in Manufacturing'
)
for(i in seq(1,10,1)){
predictions_middleclass_poly <- results_list[[i]]
lines(mobility_data_nomissings[order(mobility_data_nomissings$frac_middleclass), ]$frac_middleclass,
predictions_middleclass_poly,
col='red')
}